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Q" ■ Abstract 

The purpose of this paper is twofold: first we want to introduce a new type of hadronic interac- 
tion model (neXus), which has a much more solid theoretical basis as, for example, presently used 
models like QGSJET and VENUS, and ensures therefore a much more reliable extrapolation towards 
high energies. Secondly, we want to promote an extensive air shower (EAS) calculation scheme, 
Oh. based on cascade equations rather than explicit Monte Carlo simulations, which is very accurate 

in calculations of main EAS characteristics and extremely fast concerning computing time. We 
employ the NeXus model to provide the necessary data on particle production in hadron-air col- 
lisions and present the average EAS characteristics for energies lO^** — 10^^ eV. The experimental 
k> i data of the CASA-blanka group are analyzed in the framework of the new model. 



1 Introduction 

Although cosmic rays have been studied for many decades, there remain still many open questions, 
in particular concerning the high energy cosmic rays above 10^'' eV. One does neither know their 
composition nor the sources and acceleration mechanisms, partly because of the fact that at these 
high energies direct measurements are impossible due to the weak fiux. But since cosmic ray 
particles initiate cascades of secondaries in the atmosphere, the so-called extensive air showers, 
one may reconstruct cosmic ray primaries by measuring shower characteristics. This reconstruction 
requires, however, reliable model predictions for the simulation of EAS initiated by either protons 
or nuclei from helium to iron. The problem is that the energy of the cosmic rays may exceed by 
far the energy range accessible by modern colliders, where at most equivalent fixed target energies 
of roughly 10^^ eV can be reached. The projects of new generation EAS arrays are aimed even 
to the energy region 10^° -^ 10^^ eV, and so the gap between the existing energy limit of collider 
data and the demands of cosmic ray experiments is considerable. Moreover, the real reliable data 
limits are essentially less than above mentioned, because at present time collider experiments do 
not register particles going into the extreme forward direction and a number of other drawbacks 
may be listed. Therefore there exists a real need of "reasonable" models, implementing the correct 
physics, in order to be able to make extrapolations towards extremely high energies. 



Concerning models one has to distinguish between EAS models and hadronic interaction models. 
The latter ones like VENUS ^, QGSJET J2[ ||], and SYBILL |Q] are modeling hadron-hadron, hadron- 
nucleus and nucleus- nucleus collisions at high energies, being more or less sophisticated concerning 
the theoretical input, and relying in any case strongly on data from accelerator experiments. The 
EAS models like CORSIKA ||^ are actually simulating the full cascade of secondaries, using one of 
the above-mentioned hadronic models for the hadronic interactions and treating the well known 
electro-magnetic part of the shower. It turns out that the model predictions of EAS simulations 
depend substantially on the choice of the hadronic interaction model. In CORSIKA, the average 
electron number in EAS at primary energy 10^^ eV varies from 1.11 • 10^ to 1.62 • lO'^ (at sea 
level) depending on the hadronic interaction model |g| . So the right choice of the model and its 
parameters is extremely important. 

Recently a new hadronic interaction model neXus fj has been proposed. It is characterized 
by a consistent treatment for calculating cross sections and particle production, considering en- 
ergy conservation strictly in both cases (which is not the case in all above-mentioned models!). 
In addition, one introduces hard processes in a natural way, avoiding any unphysical dependence 
on artificial cutoff parameters. A single set of parameters is sufficient to fit many basic spec- 
tra in proton-proton and lepton-nucleon scattering, as well as in electron-positron annihilation. 
Briefly: concerning theoretical consistency, neXus is considerably superior to the presently used 
approaches, and allows a much safer extrapolation to very high energies. 

This new approach cures some of the main deficiencies of two of the standard procedures 
currently used: the Gribov-Regge theory and the eikonalized parton model, the theoretical basis of 
the above-mentioned interaction models. There, cross section calculations and particle production 
cannot be treated in a consistent way using a common formalism. In particular, energy conservation 
is taken care of in case of particle production, but not concerning cross section calculations. In 
addition, hard contributions depend crucially on some cutoff and diverge for the cutoff being zero. 

Having a reliable hadronic interaction model, one may now proceed to do air shower calculations. 
It may appear that an ideal solution from the user's point of view is to use direct Monte Carlo 
technique, where the cascade is traced from the initial energy to the threshold one and the threshold 
energy corresponds to the minimum energy registered by the array in question. But such an 
approach takes unreasonably much computer time and sometimes gives no possibility to analyze 
experimental data in the appropriate way. Already at 10^** eV to 3 • 10^® eV, the experimental 
statistics from the KASCADE experiment ||^ is ten times bigger than the corresponding number 
of simulated events. It should be mentioned that at higher energies (above 10^^ eV) it is hardly 
possible to use direct Monte Carlo calculations. Of course, some modifications, like the thinning 
method ||g], are possible and average values may be computed. 

In this paper, we follow an alternative approach, where one treats the air shower development 
in terms of cascade equations. Here, the cascade evolution is characterized by differential energy 
spectra of hadrons, which one obtains by solving a system of integro-differential equations. Crucial 
input for these equations is the inclusive spectra of hadrons produced in hadron-air collisions. These 
spectra are obtained by performing neXus simulations and finally parameterizing the results. 
Having solved the cascade equations, we finally analyse the results of the CASA-BLANKA group as 
a first application of our approach. 

2 The NEXUS Model 

The most sophisticated approach to high energy hadronic interactions is the so-called Gribov- 
Regge theory | |lO| , |ll|. This is an effective field theory, which allows multiple interactions to 
happen "in parallel", with phenomenological object called "Pomerons" representing elementary 
interactions [Q. Using the general rules of field theory, one may express cross sections in terms 
of a couple of parameters characterizing the Pomeron. Interference terms are crucial, they assure 
the unitarity of the theory. 



A big disadvantage is the fact that cross sections and particle production are not calculated 
consistently: the fact that energy needs to be shared between many Pomerons in case of multiple 
scattering is well taken into account when considering particle production (in particular in Monte 
Carlo applications), but not for cross sections ]lq |. 

Another problem is the fact that at high energies, one also needs a consistent approach to 
include both soft and hard processes. The latter ones are usually treated in the framework of the 
parton model, which only allows to calculate inclusive cross sections. 





Figure 1: The diagram representing a proton-nucleus collision, or more precisely a proton interacting with (for 
simplicity) two target nucleons, taking into account energy conservation. Here, the energy of the incoming proton 
is shared between all the constituents, which provide the energy for interacting with two target nucleons. 

We recently presented a completely new approach ||l4| , |l5| , 0| for hadronic interactions and the 
initial stage of nuclear collisions, which is able to solve several of the above-mentioned problems. 
We provide a rigorous treatment of the multiple scattering aspect, such that questions of energy 
conservation are clearly determined by the rules of field theory, both for cross section and particle 
production calculations. In both (!) cases, energy is properly shared between the different inter- 
actions happening in parallel, see fig. |l]for hadron-nucleus collisions. This is the most important 
and new aspect of our approach, which we consider a first necessary step to construct a consistent 
model for high energy nuclear scattering. 

The elementary diagram, shown as the thick lines in the above diagrams, is the sum of the usual 
soft Pomeron and the so-called semi-hard Pomeron, where the latter one may be obtained from 
perturbative QCD calculations (parton ladders). To some extent, our approach provides a link 
between the Gribov-Regge approach and the parton model, we call it "Parton-based Gribov-Regge 
Theory". 



3 System of Hadronic Cascade Equations 

EAS are produced as a result of the hadronic cascade development in the atmosphere. We charac- 
terize hadronic cascades by the differential spectra hn{E, X) of hadrons of type n with energy E 
at an atmospheric depth X, the latter one being the integral over the atmospheric density p along 
a straight fine trajectory (not necessary radial) from some point P to infinity. 



X 



p{x)dx, 



(1) 



measured usually in g/cm^ . The decrease of the average hadron numbers due to collisions with air 
nuclei is given as 



dX 



Xr, 



where the mean inelastic free path A„ (in units of mass/area) can be expressed via the average 
hadron-air cross-section a>^^^ and the average mass of air molecules niair '■ 

A« - ^ (3) 

The second process to be considered is particle decay. The decay rate in the particle cm. system 
is dhn/dr = —hn/ro, with tq being the particle life time. For a relativistic particle we find 

dhn _ ^n , . .X 



with the decay constant in energy units 



Bn = ^^, (5) 

acTo 



where to„ is the hadron mass and c the velocity of light. The ratio a — pjX depends only weakly 
on X and is for the following taken to be constant, which implies constant -B„. If we take the 
simple exponential barometric formula for the density p, we get 

PogcosO 
a = p > (6) 

-TO 

where po and Pq are density and pressure at (for example) sea level, g the gravitational acceleration 
of the earth, and 9 the zenith angle of the shower trajectory. In this paper, we only consider the 
case 9 = Q. 

Based on the above discussion, a system of integro-differential equations for the differential 
energy spectra /i„ of hadrons may be presented as 



dK{E,X) 
dX 



-K{E,x 

IE 
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™ Je 



Xn{E) EX 

Wmn{E',E) BmDjnn{E',E) 



\m{E') E'X 



(7) 
dE' 



The quantities Wmn{E' , E) and Dmn{E' , E) are the inclusive spectra of secondaries of type n and 
energy E which are produced in interactions (W) or decays (D) of primaries of type m and energy 
E' . The energy Smax is the maximum energy considered. If the type of primary hadron is Uo, its 
energy Eq and the cascade originates at Xo, then one should add 

K{E,X = Xo) = SnnJiE - Eo) (8) 

as the boundary condition. The more detailed consideration of the problem may be found elsewhere 
(see 111). 

It is reasonable to incorporate in the system nucleons (and anti-nucleons) , charged pions (Stt = 
114 GeV), charged kaons (Bk = 852 GeV), and neutral kaons (Sk£ = 205 GeV). As Bko = 
1.19 • 10^ GeV, there is no sense to account for neutral kaons K^ at energies <C Bk°- But these 
particles should be included if their energy exceeds O.OIBk". The values of decay constants are 
given at the height 11 km. 

The computational technique to solve the system (|^) is based on the same principle as the 
traditional approaches ||l^, |l^, but some improvements are introduced, which enable one to avoid 
too small steps when integrating over the depth ||l^. One discretizes the energy as 

(9) 



with i?min — 1 GeV and c such that the number of points per order of magnitude is 10 -^ 20. 
Replacing the integral in the right-hand side of (0) by the corresponding sum, one may write 



dX 



= -Ku{x) 



Bn 



^ £;!,„, (X) 



m j„ 






(10) 



with hni{X) =hn{Ei,X), A™ = A„(i?i), and 
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E, 



(11) 
(12) 



with i?min(J) = E^l^fc and E^aAi) ^ Ei ■ y/c { ii i ^ j then Emax{i) = -E^ ). The factor iS/i;^ 
has been added in the integral to ensure exact energy conservation. The spectra WmniEj, E) and 
Dmn{Ej,E) must be provided in order to calculate the above integrals. The standard recipe for 
calculating Dmn{Ej,E) may be found elsewhere ( see ||l^ or |g] ) and impHes no difficulties. The 
calculation of Wmn{Ej,E) must be based on neXus simulations, as discussed below in detail. 
Since the Ej dependence is very smooth, there is no point to calculate the spectra for all energies 
Ej. One rather chooses some reference energies Ej/ to calculate the spectra based on NeXus 
simulations, and then interpolates to obtain the spectra for the other energies. This method has 
proven to be superior concerning the computational time when compared with the direct spectra 
calculations over all energies. For the calculations in this paper we choose two reference points 
per order of magnitude, starting from 10^^ eV: Eji — 10^^ eV, 10^^'^ eV, 10^^ eV etc. up to 10^'^ eV 
which is at present the maximum energy attainable in the NeXus model. As the NeXus model is 
not valid at energies below 10^^ eV, the data obtained with other codes must be borrowed. In this 
work we employ results obtained in [gy|, |gl| which are close to predictions of the GHEISHA code 
^^ used in CORSIKA. 

The solution of the homogeneous equation 



dK,{X) 
dX 



-hm{X) 
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Bn 

eJc 



(13) 



has the form 



hni{Xo + AX) = hni{Xo) exp 



X, 



Xm J \Xo + AX 



B„/Ei 



if for simplicity we neglect the weak dependence of Bn on X . The solution of the full equation 
([icl) may be written as 



hm{Xo + AX) 



/.™(X„)exp^-^}(^-|^ 



B„/Ei 



(14) 



EE 



m j„ 



Xa+AX 



exp 




B„/Ei 



dX' 



which may be verified directly. The above formula allows to calculate hni at depth Xg + AX, 
provided that hni{X) at X = Xq is known and all hmj{X) are also known for Xg < X < Xq + AX 
and j > i. So, starting from hni{Xo) one may sequentially find hni{Xo + AX) and so on. Test 



calculations show that for Ei > Bn/3 it is quite sufficient to use Simpson's formula. Here, one 
needs the values for hmj{Xo + AX/2), which are obtained via interpolation. Whereas for pions 
this works without problem, for kaons the requirement Ei > Bj^±/3 becomes more restrictive and, 
what is especially essential, errors for one component manifest themselves in other components. In 
order to retain accuracy without resorting to excessively small AX, the integration over X' in eq. 
(Ill) is approximated as 



Xa+AX 



I 



fix') 



X' 



Xo + AX 



B„/Ei 



dX' = Ai f{Xo) + A2 f{Xo + AX/2) + A3 f{Xo + AX) (15) 



and coefficients Ai, A2, A3 are found from the condition that (|15|) is exact for a second order 
polynomial. The accuracy of ~ 1% may be achieved with AX ~ 5 g/cm^. 

Other EAS characteristics (e.g. electron and muon numbers) are computed in a traditional 
way as corresponding functionals from functions hn{E, X). Usually one assumes that neutral pions 
decay immediately at the generation point and do not contribute to the development of the hadronic 
cascade. This assumption is quite adequate at energies below the corresponding decay constant 
which for neutral pions is about 3 • lO^^eV. Moreover, as primary particles are nucleons there is an 
additional factor of ~10 in our favour. So the number of neutral pions produced at depth X may 
be obtained as 
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where we substitute n by tt" in the second term of the right-hand side of eq. (|ic 
number Ne at depth T is given as 



(16) 



The electron 



Ne{T) 



i=l 



Ko{E„Xa)NGiy,t) 
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■dXo 



(17) 



where y = hi{Ei/l3) is the logarithmic energy in units of the critical energy of electrons in air 
(/3 = 81.10^ eV) and t = {T - Xo)/To is the depth difference in radiation units (Tq =37.1 g/cm^). 
The factor 2/(1 + s) accounts for energy sharing between two photons, with s being the shower 
age parameter s ~ 3t/{t + 2y). The function Nq is referred to as Greisen's formula p^ . 



0.31 

No = —pr exp 



1 - - ln(s) 



(18) 



which predicts the electron number at depth T in the shower, produced by a primary photon with 
energy Ei at depth Xq. 

As a rule, experimental EAS arrays can detect muons with energies above a certain threshold 
Et\i-[p,. The number of such muons may be obtained as follows: 



N^{E^ > i?thr^, T) = J2 E E / ^™^(^) Jt ^"/^ ^(^^' ^' ^*''^'^' ^) '^^' ^^^^ 

^vi . . 17 ^ IT- „■ \ ,- -^ J 



m i:Ei>Eth,^, j>i 



where D^^^ defines the number of muons with energy Ei resulting from the decay of hadron m with 
energy Ej, W{Ei,X,Ethiii,T) is the probability that a muon produced at depth X with energy 
Ei will survive between X and T and its final energy at T will be greater than i^thr^. 



For the most important channels of muon production (tt^ 



A* 



and if ^ 



responsible for about 95% of all muons) values of D-j^ are governed by the simple two body decay 



kinematics. If we consider the atmosphere to be isotermic then the function W{Ei,X,Ethr^l,T) 
may be written explicitly: 

WiE,,X,E,^,„T) ^ e (e., - g,,.., - (T - X)gj || ^' ^'^^ ^^'^^ I (20) 

where 9 is the step function and ^ is the ionization loss rate. 

4 Calculations of the Inclusive Spectra 

In this section, we discuss how to obtain the matrices W^^^ representing inclusive particle spectra, 
based on NeXus calculations. 

The first step evidently amounts to performing a number of NeXus simulations to calculate 
the energy spectrum of produced (secondary) hadrons of type h, 

E^^{E,^,E), (21) 

in a hadron plus nitrogen reaction /iinN — > hX for different incident hadrons hin, each one at 
different energies Ei^. For the incident energies, we use 

E(^ = 10^+1 /^GeY, J =0,1,... 12, (22) 

and for the incident hadron types, we take nucleons, charged pions, charged kaons, and neutral 
kaons Kl , for the secondary ones in addition neutral pions and Kg (spectra of Kg are assumed to 
be identical to Kl). All other hadrons are assumed to decay immediately. 

neXus is based on the Monte Carlo technique, implying automatically statistical fluctuations, 
which are in particular large, when the production of secondaries with energies approaching the 
primary one is examined. The influence of the limited statistics obtained via NeXus Monte Carlo 
simulations may be eliminated to some extent if an appropriate smoothing procedure is applied. 
Such a procedure provides the opportunity to exploit the corresponding 'smoothened' spectrum 
at a flxed primary energy as a continuous function of the energy E of secondaries, instead of 
considering discrete values only. Moreover, it enables one to impose certain restrictions on the 
shape of inclusive spectra. For example, we assume the E dependence of W^ for a: = E/E\n -^ 1 
to be proportional to a;~^(l — x)", where a is taken from theoretical considerations or just as a flt 
parameter. We use the Levenberg-Marquardt (LM) method |g4| to flt the Monte Carlo spectra by 
analytic continuous functions in E, 

E^^^{El,E) ^ W^,,{El,E) (23) 

for the different incident hadrons hm and secondaries /i„ for the above-mentioned values El^ for the 
incident energy. This method flnds the minimum y^ and its algorithm consists in the combination 
of the inverse-Hessian method and the steepest descent method. It is one of the standard non-linear 
least-square approaches and its detailed description may be found elsewhere ||2^. The statistics 
of 10^ events for each of the reference energies El^ is sufflcient to calculate Wmn{El^,E) within 
~ (0.1 -=- 0.2)% accuracy. 

As an example, we show in flgs. g and || the Monte Carlo results for incident charged pions 
with energies of 10^^ eV and 10^^ eV, together with the corresponding flt functions Wmn{El^,E). 
Similarly excellent flts are obtained for all the other spectra. For purposes of more efflcient com- 
puting, two sets of Monte Carlo spectra are used, which exploit different £'-scales: a linear one for 
large energies of secondaries (right flgures) and a logarithmic scale for small ones (left flgures) . 
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Figure 2: Inclusive spectra (points) from neXus simulations and the corresponding smoothened spectra 
WmniE?^, E) (lines) for incident charged pions and different secondaries h„, for incident energy Ein = 10^^ eV, 
as a function of the secondary energy E (left) and x = E/E-^^ (right). 
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Figure 3: Same as fig. kI but for incident energy Bin = 10^^ eV. 



Using analytic expressions for the spectra, we can now proceed to calculate integrated spectra, 
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W^n{El.E)dE 



with Ei — E'min c' for 1 < i < in 



For the actual calculation, we use c = 10" and i 



(24) 



81. 



We obtain the integrated spectra for arbitrary incident energies i?in via an interpolation formula 
W'4„(£'in), and we can thus calculate these quantities in particular for all the energies Ej-. 






(25) 



One could have calculated the integrated spectra directly for all the energies Ej^ but due to the 
weak energy dependence of the spectra it is much more efficient to proceed as discussed above. 

5 Solving the Cascade Equations 

Having all the ingredients, we are now able to solve the cascade equations as discussed above to 
obtain the differential hadron spectra hn{Ei,X). The method adopted gives the possibility to 
calculate average characteristics of EAS within ~ 1% accuracy. Based on the inclusive spectra, we 
may calculate numbers of different hadrons TV/j. as well as numbers of electrons N^ and Nf^ for a 
number of observation levels X in the atmosphere. In fig. Q, we show the dependencies of the 
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Figure 4: Number of electrons A/^e, muons Af^ (E^ > 1 GeV), and all hadrons A^j, (E^ > 50 GeV) as a function of 
the atmospheric depth X for different incident energies -Ein (in eV, from top to bottom): 10^®, 10^ 



lO^'^ 10^ 



number of electrons Ne , muons 7V^ (i5^ > 1 GeV) , and all hadrons Nh {Eh > 50 GeV) on the depth 
X for different incident energies. One observes the expected increase of the particle numbers with 
energy and as well the shift of the shower maximum towards larger X with increasing energy. In 
fig. 0, we show the corresponding characteristics for individual hadrons. The pions are by far most 
dominant, followed by nucleons, then charged kaons and K^, whereas Kg are the least frequent due 
to the short life time. In fig. ^, we show the total hadronic energy as a function of the atmospheric 
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Figure 5: Numbers of different hadrons A''^. (B^. > 50 GeV) as a function of the atmospheric depth X for different 
incident energies E\^. 
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depth X for different incident energies. Obviously the hadronic energy is highest for the highest 
incident energy. With increasing atmospheric depth, the hadronic energy drops exponentially, due 
to its conversion into the energy of electro-magnetic cascade (and to some extent into muon and 
neutrino energy). 

6 Some Results of Calculations 

It is certainly the main purpose of this paper to introduce a new, sophisticated hadronic interaction 
model, particularly suited for high energy hadronic interactions, and, at the same time, to explain 
the use of cascade equations rather than explicit simulations for the air shower calculation. So we 
do not want to present extensive appHcations of this approach, but rather discuss one instructive 
example. 

The CASA-BLANCA group Ig^ published recently very interesting results concerning the com- 
position of cosmic rays in the energy range 10^''-10^^ eV. From the results shown in the previous 
section, we can easily calculate the shower maximum Xmax as a function of the incident energy 
Ein- We calculate as well the shower maximum for incident nuclei (with mass number A) assuming 
that it is given by the result for nucleon at a reduced energy Ein/A. As it had been shown (see, 
for example, ||2|, ^) this simple superposition prescription works well for average characteristics of 
nucleus-induced EAS. In fig. we show the shower maximum X^aa^ for incident nucleons (upper 
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Figure 7: The shower maximum Xmax for incid 
the incident energy Ei^ together with the data || 



at nucleons (upper curve) and iron (lower curve) as a function of 



curve) and iron (lower curve) as a function of the incident energy Em together with the data. 
Comparing our calculations with the data we find an excellent agreement with the results of anal- 
ysis carried out in [Q on the basis of QGSJET and VENUS models, i.e. we confirm the change to 
a heavier primary composition at energies above the knee of the primary cosmic ray spectrum, 
which has been earlier reported by the group of Moscow State University |29, BOl, as well as the 
novel feature discovered by the CASA-BLANCA group - "lightening" of the composition at energies 
just before the knee. 

The above statement can be made more quantitative by studying the so-called mean nuclear 
mass {A), which is defined to be the nuclear mass which would fit best the experimental data. In 
fig. ^, we plot the mean logarithmic mass \n{A) as a function of the incident energy E'in together 
with QGSJET and VENUS results obtained in ||2^. As already mentioned above, the mass number 
has clearly a minimum at around 10^^-^ eV, for higher energies the mass is increasing again. The 
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neXus results are quite similar to the QGSJET ones, whereas VENUS has a tendency towards higher 
masses. One should keep in mind, however, that VENUS is strictly speaking already outside its 
energy range of validity, which makes its prediction somewhat uncertain. 

It is worth to note that the observed qualitative behaviour for the primary composition may 
be obtained in the framework of the diffusion model for cosmic ray propagation if one assumes a 
large magnetic halo for the Galaxy with antisymmetric magnetic field |^. As it has been shown 
in ll^ such a configuration of the galactic magnetic field agrees with the measurements and allows 
to explain naturally the observed "sharpness" of the knee in the primary energy spectrum. 
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Figure 8: The mean logarithmic mass In(^) as a function of the incident energy Ei^. We compare our neXus 
results with the qgsjet and venus ones obtained in |26|. 



7 Summary and Outlook 

We introduced a new type of hadronic interaction model (neXus), which has a much more solid 
theoretical basis when compared, for example, to presently used models such as QGSJET and VENUS, 
and so provides much more reliable predictions at superhigh energies where there are no collider 
data yet. A particular feature of the model is a consistent treatment for calculating cross sections 
and particle production considering energy conservation strictly in both cases. In addition, one 
introduces hard processes avoiding any unnatural dependence on artificial cutoff parameters. Using 
a single set of parameters, NeXus is able to fit many basic spectra in proton-proton and lepton- 
nucleon scattering, as well as in electron-positron annihilation. So it is worth to point out once 
more that concerning theoretical consistency, NeXus is considerably superior to the models in 
current use, and allows a much safer extrapolation to superhigh energies which are very important 
in cosmic ray studies. 

We explained in detail an air shower calculation scheme, based on cascade equations, which is 
quite accurate in calculating main characteristics of air showers, being extremely fast concerning 
computing time. We perform the calculation using the NeXus model to provide the necessary 
tables concerning particle production in hadron-air coUisions. 

As an application, we calculated the shower maximum as a function of the incident energy for 
incident protons and iron, to compare with corresponding data. Based on these data, we calculated 
as well the so-called mean nuclear mass. We are thus able to confirm that the average mass as a 
function of the incident energy shows a minimum around 10^^'^ eV. 

Problems where it is possible to ignore fluctuations are not too numerous. So an unavoidable 
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question arises how to retain calculation efficiency at reasonable level and at the same time to 
account for fluctuations. But the answer to this really crucial question is well known and proves 
to be rather simple. One needs to employ a combination of the Monte Carlo techhnique and 
numerical solutions of hadronic cascade equations in the atmosphere. The explicit simulation 
should be carried out from the primary (initial) energy Em to some threshold Sthr = kEo, where 
k ~ 10~^ — 10~^. Numerous calculations showed (see, for example, [Q) that there is no sense in 
going below this threshold as it does not increase the accuracy of EAS fluctuation determination. So 
contributions from hadrons with energies below Ethi may be accounted for in average. Numerical 
solutions of hadronic cascade equations can produce corresponding tables with sufficient accuracy 
and in a very short time. 

Here we do not consider calculations of lateral distributions of shower particles as this problem 
will be discussed in our next paper. But we would Hke to point out that there is a rather important 
class of problems connected with giant EAS arrays aimed at primary energies > lO^^eV (see fs^ ). 
For these arrays it is sufficient to calculate lateral distributions only at large distances (above 100 
m) from the shower axis. Such a situation makes it possible to treat the problem as a combination of 
the one-dimensional approach for hadrons with energy > lO^^eV and the rigorous three-dimensional 
technique for low energy region only. 
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